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The linear response of an isolated, homogeneous granular fluid to small spatial perturbations is 
studied by methods of non-equilibrium statistical mechanics. The long wavelength linear hydrody- 
namic equations are obtained, with formally exact expressions for the susceptibilities and transport 
coefficients. The latter are given in equivalent Einstein-Helfand and Green-Kubo forms. The context 
of these results and their contrast with corresponding results for normal fluids are discussed. 
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Granular fluids consist of mesoscopic particles acti- 
vated by external forces or initial conditions, and appear 
ubiquitously both in nature and in technological appli- 
cations. Their flow properties are remarkably similar to 
those of normal fluids in many cases [l]. This has led 
to the speculation that the methods of non-equilibrium 
statistical mechanics developed for normal fluids might 
be adapted to explore the properties of granular fluids 
[^1^. To date, this has been realized in the examples 
of numerical molecular dynamics simulation and kinetic 
theory. However, the extensive formal theoretical tools 
have been applied to a more limited degree (e.g., impu- 
rity diffusion 4] and liquid state structure [5j)- In some 
respects experimental studies of the response functions of 
interest here are more advanced than those of theory (see, 
for example, [^). Here we report an application of linear 
response methods to obtain the hydrodynamic descrip- 
tion of small perturbations in an isolated, homogeneous 
granular fluid. Although this example is perhaps of lim- 
ited direct experimental relevance, it provides an instruc- 
tive controlled test for the extension of these methods to 
granular fluids, illuminating important differences from 
normal fluids. 

The simplest model of a granular fluid is considered: N 
smooth, inelastic, hard spheres (d = 3) or disks [d — 2) of 
mass m. The only difference from a corresponding model 
for normal fluids is a loss of energy in each binary colli- 
sion, characterized by a restitution coefficient < a < 1, 
with a = 1 corresponding to the elastic limit. At the mi- 
croscopic level, the fluid state is specified by the positions 
and velocities of all particles, i.e. by a point in the 2dN- 
dimensional phase space T = {qi, - ■ ■ ,qj^,vi, - ■ ■ ,vj^}. 
This point evolves according to a deterministic dynamics 
of free streaming and collisions. The initial macroscopic 
state is specified by a probability density p(r,0), char- 
acterizing the known (incomplete) initial data. Conse- 
quently, all the ingredients required for a statistical me- 
chanical description are available. In particular, the Li- 
ouville equation determines the macrostate p (T, t) at all 
later times. The elements of a hydrodynamic description 



then follow from averaging the local number, energy, and 
momentum densities to obtain the exact macroscopic bal- 
ance equations. From them, the phenomenological hy- 
drodynamic equations are obtained with the additional 
assumption of constitutive equations for the fluxes, as 
well as for the source term appearing in the equation 
for the energy. For systems with small gradients of the 
hydrodynamic flelds, the Navier-Stokes equations for a 
granular fluid are obtained. Of course, the transport co- 
efficients in these equations are unknown and must be 
supplied separately from experiment. 

For a normal fluid, a formal theoretical description of 
hydrodynamics has been given by considering the linear 
response of an equilibrium fluid to small spatial pertur- 
bations 0] ■ The solutions to the linearized hydrodynamic 
equations, the so-called hydrodynamic modes, can then 
be identified from the long time, long wavelength be- 
havior of the response functions. As a result, formally 
exact expressions for the transport coefficients are ob- 
tained in terms of time correlation functions. These are 
the Einstein-Helfand (EH) and Green-Kubo (GK) forms, 
which have proved useful as the basis for practical (ap- 
proximate) theoretical predictions of these coefficients. 
The objective here is to report the application of this 
linear response approach to granular fluids. 

The isolated granular fluid does not support a sta- 
tionary, homogeneous equilibrium state. Instead, the 
simplest homogeneous state, known as the homogeneous 
cooling state (HCS), changes in time due to collisional 
energy loss. Spatial perturbations of the HCS therefore 
have a dynamics due to both the perturbation and the 
background reference state. However, the latter occurs 
only through the average energy and can be removed by 
a suitable choice of dimensionless variables. Effectively 
then, the analysis is similar to that for the equilibrium 
case, although the reference state is not the Gibbs distri- 
bution. The distribution function of the HCS, po{T), is 
an example of a "normal distribution" , in the sense that 
its time dependence occurs entirely through the temper- 
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ature, 

Po{T,t) = {iv^it))-^'' pl {{q^J/^,v,/vo{t)},na) . (1) 

Here = qi qj and £ is the mean free path. The uni- 
form hydrodynamic fields are the average number den- 
sity, no, and the temperature, To{t), the latter occurring 
through the thermal velocity VQ{t) = ^2TQ{t)/m . The 
granular temperature is defined in terms of the average 
energy density, cq, in the same way as for equilibrium, 
Co — (inoTo/2. The HCS distribution is the solution to 
the Liouville equation in the form 



Cpo = 0, £ = L ' ^° 
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where L is the Liouville operator for hard spheres |^|^, 
and Co = — c^tlnTo is the "cooling rate". 

To study the response of the HCS to spatial pertur- 
bations, an initial local HCS, p;(r,0), is chosen. This 
distribution is a functional of specified nonuniform hy- 
drodynamic fields n(r), T{r), and u{r), characteriz- 
ing the average number density, energy density, and 
momentum density for this ensemble. To construct 
it, consider first the HCS for an inhomogeneous fluid 
(i.e., the solution to ^ in an external inhomogeneous 
field). Since the average nonuniform density is a func- 
tional of the external field, inverting that relationship 
gives the HCS distribution as a functional of the density 
pi^nh) _ (^gy^'j-d.N p* (^!^q^-/£^i,i/vo}\n). The special case 
in is the value of this functional at the uniform state. 
Next, the non-uniform flow field and temperature are in- 
corporated through the replacements of Vi by Vi — u{qi) 
and Vo{Tq) by Wo(Q'i) = ''Jo[T{q^i)] for each particle in the 
arguments of pg™''' . The local HCS distribution is then 

(3) 

This initial preparation is natural, in the sense that it 
is expected that locally there is a rapid relaxation of the 
velocity distribution as each region attempts to approach 
the HCS at its own values of the hydrodynamic fields, ft 
is the analogue of rapid coUisional approach to a local 
equilibrium state for normal fluids. For smooth spatial 
variations of small amplitude, to linear order 



p(0) = pi{0)~po + Y, J dr 
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where the y;3(r, 0) are linear combinations of the hydro- 
dynamic fields, as specified below. The second equal- 
ity follows by considering a single Fourier component for 
these fields, and it defines 4)p{k). 



This preparation is also special in the sense that it ex- 
cites only hydrodynamic modes in the long wavelength 
limit. This follows from the fact that the yp can be cho- 
sen such that = (j)p{k = 0) are eigenfunctions of the 
extended Liouville operator C 
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with the eigenvalues {\p} = {0, Co/2, -Co/2, .., -Co/2}. 
The specific set of fields in Q leading to this result are 

+2|;^i^^,^,^|. (6) 

no To a In no no vq Vq J 

The field fiow u has been decomposed into its longi- 
tudinal component in the direction of fc = fc/fc, and 
the remaining d — 1 transversal components denoted by 
u±. The eigenvalue —Co/2 is d-fold degenerated. Sig- 
nificantly, the eigenvalues A/j are the same as those of 
the linearized macroscopic balance equations for these 
variables in the long wavelength limit. Hence they can 
be considered as hydrodynamic modes reflected in the 
microscopic dynamics. These observations are critical 
to the following linear response analysis. In the elastic 
limit, the eigenfunctions become linear combinations of 
the global conserved number, energy, and momentum. 

The response of the hydrodynamic fields Syp to the 
initial perturbation Q is 



{k* 



Cp^ik' 
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Cp^{k\s)5y*{k* ds 



dt. 



,s)^ J dr*SY*{k*)e-^^'r^{-k). 



(7) 
(8) 



The dimensionless time s is a measure of the average 
number of collisions, and Yp is the phase function whose 
average is yp. Moreover, the asterisks indicate dimen- 
sionless variables in the units defined by £, vo, and m 0]. 
Clearly, if the 5y*p{k* ^ s) obey hydrodynamic equations at 
long space and time scales, the corresponding excitations 
must exist in the spectrum of the correlation functions 
Cpj{k*,s). This suggests defining an exact transport 
matrix JC^pik* ts) by 

[d,+IC{k*,s)]Sy*{k*,s)^0, (9) 

with the identification 

IC{k*,s) = - [dsC{k*,s)]C-\k*,s). (10) 

A matrix notation has been introduced for simplicity. 
The hydrodynamic equations are then identified from the 
matrix /C(fe*,s) for large s and small k*. In particular, 
the linear Navier-Stokes approximation is defined by an 
expansion of IC{k*,s) to order A:*^, with the coefficients 
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evaluated for large s. In general, this procedure does not 
lead to tractable results since initial perturbations excite 
microscopic as well as hydrodynamic modes. However, 
the special initial preparation chosen here is such that 
only hydrodynamic excitations occur in the long wave- 
length limit. This follows directly from jSJ and © 



Pi- 



(11) 



As a consequence, the expansion of H10|l is assured to 
start on the hydrodynamic branch. Define 

C (fc*, s) = e-'^ + ik*C^^^ (s) + {ik*f (s) + • • • , 

(12) 

where A^^ — \*^5afi- Then it is found 

lC{k*, s) = A + ik*K:^^\s) + {ik*fK^'^\s) + ■ • • , (13) 
with 



/C«(s) = -[(9,+A)C(i)(s) 



(14) 



/C(2) (s) = - \{d, + A)C(2) (s) + /C'l^ (s)C(i) (s) 



(15) 

The transport coefhcients and susceptibilities of the phe- 
nomenological Navier-Stokes equations can be identified 
with appropriate matrix elements in these expressions: 
gives the Euler coefficients, A^*-^-* gives the Navier- 
Stokes coefficients, and so on. In particular, these ex- 
pressions provide an extension of the EH representation 
of transport coefficients to granular fluids. The GK forms 
follow from these results by partial integration. The de- 
tailed expressions and analysis will be given elsewhere. 

In this brief presentation, further details will be lim- 
ited to the case of the transverse flow field response, i.e. 
6y*^{k*,s) associated with one component of in 
The only non- vanishing element Cai_ is Cx^. Also note 
that C^ll = fr om symmetry. Then, Eq. ^ becomes to 
order fc^ 



5yl{k\s) = Q. (16) 



If /C^j^(s) has a limit for large s, this becomes the shear 

diffusion equation, with 1C\_\_ identified in terms of the 
dimensionless shear viscosity rj* = r]/mnQ£vQ, 



lim IC^liis) 



S>>1 



= lim ds — 

s>>i N 



lim ds 

S>>1 



'in 



dV* 



M* 



U*pl, (17) 



where —C* is the adjoint of C . The phase function M* 
is the space moment of the momentum density 



M* 
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where ei is the unit vector in the direction of u±_ consid- 



ered. The function 7V*pj5 



a(i) 



is given by 



N 



dr* [k ■ r*] ei 
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(19) 



In the elastic limit and for po Pe, the equilibrium 
canonical distribution function, the Helfand form for the 
shear viscosity of a normal fluid is recovered 

Vli^s, = - hm ds^J dV* (e^^'M*) M*p:. (20) 

Both (|17|l and H2U|) contain correlation functions involv- 
ing the momentum density moment M. For a normal 
fluid, this is an autocorrelation function, while for the 
granular fluid a new moment appears that is deter- 
mined from the HCS distribution. 

These correlation functions must grow as s for long 
times to obtain a constant viscosity. This is analogous 
to the definition of the diffusion coefficient for a Brow- 
nian particle in terms of its mean square displacement, 
and for that reason the representation above is referred 
to as the EH form. An equivalent representation is ob- 
tained by carrying out the time derivative, shifting the 
time dependence, and integrating with respect to time. 
The result is 



V 



1 J dr* F*N*pl 
lim [ ds'— [ dT*F 



s<l 



N 



V ' F Or. 



The fluxes are given by 



(22) 

This is the GK representation. For normal fluids, F* 
becomes the volume integrated momentum flux and the 
second term of (|21|l becomes the expected time integral 
of a flux autocorrelation function. As noted below, hard 
sphere systems entail the additional first term on the 
right side even for normal fluids. Moreover, as with the 
EH form, the correlation functions for the granular fluid 
involve new phase functions defined from the HCS distri- 
bution. 

The objective here has been to explore the application 
of linear response methods to the simplest model and 
state of a granular fluid. The details for all hydrodynamic 
modes and all transport coefficients will be reported else- 
where. These results illustrate the applicability of such 
methods and their potential for more complex conditions 
of experimental interest. They also expose many sub- 
tleties in the differences from corresponding results for 
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normal fluids and provide a warning against a simple 
translation of normal fluid theory to granular fluids. In 
closing, some observations are highlighted: 

i) The identification of the transport matrices in H13|l - 
()I5|1 constitutes a formal derivation of the linear Navier- 
Stokes equations. Important limitations on this deriva- 
tion are placed by restrictions on length and time scales. 
As for normal fluids, it is expected that the conditions 
are times long compared to the mean free time (s ^ 1) 
and wavelengths long compared to the mean free path 
{k* = H < 1). 

ii) These formal results have no a priori limitation on 
the density or degree of inelasticity. However, peculiar- 
ities of the isolated system (cluster instability [lO], in- 
elastic collapse must be confronted at high densities 
and large inelasticity. These problems can be removed in 
practice by considering small system sizes and protocols 
for simulations to avoid dangerous configurations with no 
physical consequences. 

iii) The transport coefficients are given in terms of av- 
erages over the HCS reference state, which differs from 
the equilibrium canonical distribution. Consequently, 
one of the two phase functions in the EH and GK ex- 
pressions is different from that for normal fluids. 

iv) The generator L for hard sphere dynamics is singu- 
lar and must be expressed in terms of binary collision op- 
erators rather than forces. Important structural changes 
result in the form of the GK expressions and the asso- 
ciated microscopic fluxes The GK expressions for 
both elastic and inelastic hard sphere fluids have a new 
contribution due to instantaneous collisions, as well as 
the usual time integral of correlation functions of fluxes. 
The latter are expressed in terms of the binary collision 
operators and are different for forward and backward dy- 
namics, both occurring in the GK expressions. 

v) The reference HCS is time dependent, due to col- 
lisional cooling. A dimensionless form of the equations 
is required to eliminate this dependence, resulting in a 
new generator C rather than the Liouville operator L. 
In this representation, the HCS is a stationary solution, 

vi) The stationary HCS distribution is singular in the 
sense that it is restricted to a zero total momentum P 
surface, i.e., /Oq (x 5{P*). Furthermore, the elastic limit 
is proportional to a delta function in the total energy, 
since an isolated system is being considered. This con- 
stant momentum, constant energy ensemble is often re- 
ferred to as the "computer" ensemble, as it corresponds 
to usual conditions in molecular dynamics simulation. 
For granular fluids, the energy changes, even in isolation, 
but the constant momentum condition must be retained 
for consistency with C = Q. 

vii) The presence of hydrodynamic excitations is ex- 
pected for a wide class of initial perturbations. However, 
the choice made here is both physically and mathemat- 
ically motivated. As the HCS distribution characterizes 



the homogeneous state through its dependence on the 
hydrodynamic fields, it is natural to consider their spa- 
tial perturbation through this same dependence. More 
important for the analysis here is that the resulting per- 
turbation excites only hydrodynamic modes in the long 
wavelength limit. This follows from the property jSJ), 
and plays a central role in deriving the explicit forms for 
the expansion H13(l . For general initial preparations, the 
leading term in will not be A, but rather some com- 
plex time dependent quantity that must be calculated to 
be A at long times. 

viii) An earlier derivation of GK expressions leads to 
different results The method is formally correct, 

but does not have the simplifying features of the initial 
preparation considered here. In addition, some of the 
special properties of the hard sphere generator were not 
accounted for in detail, so that different fluxes are ob- 
tained and the instantaneous contribution is missed. 
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